Contents

1 Packagees

library(Vennerable)
library(UpSetR)
library(FDb.InfiniumMethylation.hg19)
library(tidyverse)
library(reshape2)
library(limma)
library(GenomicRanges)
library(DESeq2)
library(pROC)
library(ggplot2)
library(beeswarm)
library(plyr)
library(jyluMisc)

2 Functions

labeler <- function(x){
  x[x==0] <- "Res.WT"
  x[x==1] <- "Sens.WT"
  x[x==2] <- "Mut"
  factor(x)
}

diff.cont <- function(df, cnt){
  lev <- levels(factor(cnt))
  design <- model.matrix(~0+factor(cnt))
  colnames(design) <- c(as.character(lev[1]),as.character(lev[2]))
  contrast <- makeContrasts(g0 - g1, levels = design)
  fit <- lmFit(t(df), design)
  fit <- contrasts.fit(fit, contrast)
  fit <- eBayes(fit, 0.01)
  tT <- topTable(fit, adjust="fdr", sort.by="B", number=10000)
  tT
}

diff_count <- function(df1, df2, f){
  df1 <- df1[rownames(df1) %in% rownames(df2),]
  f_temp <- df2[rownames(df1), f]
  tt <- diff.cont(df1, paste0("g",f_temp))
}

top_sel <- function(tt, frac = 0.1){
  rownames(tt)[head(order(tt$adj.P.Val), (dim(tt)[1]*frac))]
}

enricher <- function(DEres, name, gmt){
  to.ret <- list()
  res <- data.matrix(DEres)
  res <- data.frame(res)
  res$symbol <- rownames(DEres)
  res <- res %>% filter(P.Value < 0.05)
  res$P.Value <- res$P.Value * sign(res$logFC)
  res <- res %>% dplyr::select(P.Value, symbol) 
  colnames(res)[1] <- "stat"
  res <- res %>% arrange(desc(stat))
  rownames(res) <- res$symbol
  res$symbol <- NULL
  rg <- runGSEA(inputTab = res, gmtFile = gmts[[gmt]])
  # rg$Name <- gsub(paste0(name, "_"),"", rg$Name)
  rg <- list(rg)
  names(rg) <- name
  to.ret$rg <- rg
  to.ret$enBar <- plotEnrichmentBar(resTab = rg, pCut = .2, ifFDR = FALSE, setName = gmt)
  to.ret
}

path.sel <- function(res){
  colnames(res)[c(4,6)] <- c("p.up", "p.down")
  res %>% dplyr::filter(p.up < 0.2 | p.down < 0.2) -> res
  res$Name
}

3 Load Preprocessed Data

load("../data/metadata/patmeta_180504.RData")
load("../data/Objects/exprMat.RData")
load("../data/Objects/CPS1000.RData")
load("../data/Objects/ic50.RData")
load("../data/Objects/meth.RData")

4 Datasets

vennList <- Venn(list(metadata = patMeta$Patient.ID[!is.na(patMeta$TP53)], RNAseq = colnames(exprMat), CPS1000 = rownames(cps.data), methylation = rownames(meth.df)))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))

upset(fromList(list(metadata = patMeta$Patient.ID[!is.na(patMeta$TP53)], RNAseq = colnames(exprMat), CPS1000 = rownames(cps.data), methylation = rownames(meth.df), IC50 = rownames(ic50.data))), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "orange")

5 Preprocess Data

patMeta <- data.frame(patMeta)
patMeta <- patMeta[!is.na(patMeta$TP53),]
rownames(patMeta) <- patMeta$Patient.ID
colnames(patMeta) <- make.names(colnames(patMeta))
colnames(cps.data) <- make.names(colnames(cps.data))
drug_meta <- merge(patMeta, cps.data, by = 0)
rownames(drug_meta) <- drug_meta$Patient.ID
exprMat <- data.frame(exprMat)

6 Define Cut-off for Nutlin Resistency

my_roc <- roc(predictor = drug_meta$Nutlin.3a, response = drug_meta$TP53)
cut.off <- coords(my_roc, "best", ret = "threshold")
cut.off
## [1] 0.8670314
plot.roc(
  predictor = drug_meta$Nutlin.3a,
  x = drug_meta$TP53,
  print.auc = TRUE,
  auc.polygon = F,
  grid = c(0.1, 0.2),
  grid.col = c("green", "red"),
  max.auc.polygon = TRUE,
  auc.polygon.col = "blue",
  print.thres = TRUE)

outliers <- rownames(drug_meta)[which(drug_meta$TP53 == 0 & drug_meta$Nutlin.3a >= cut.off)]
  
drug_meta$nut <- ifelse(rownames(drug_meta) %in% outliers, 1, 0)

beeswarm <- beeswarm(
  Nutlin.3a ~ TP53,
  data = drug_meta,
  method = 'swarm',
  pwcol = as.numeric(drug_meta$TP53) - as.numeric(drug_meta$nut),
  do.plot = F
)[, c(1, 2, 4)]
colnames(beeswarm) <- c("x", "y", "group")
beeswarm$group <- labeler(beeswarm$group)


beeswarm.plot <- ggplot(beeswarm, aes(x, y)) +
  xlab("") +
  scale_y_continuous(expression("Nutlin.3a Averaged")) + 
  geom_boxplot(aes(x, y, group = ifelse(group %in% c("Res.WT", "Sens.WT"), "WT", "Mut")), outlier.shape = NA) +
  geom_point(aes(colour = factor(group))) +
  scale_x_continuous(
  breaks = c(1:2),
  labels = c("unmutated", "mutated"),
  expand = c(0, 0.3)
  ) + geom_segment(aes(
    x = 0.8,
    y = cut.off,
    xend = 1.2,
    yend = cut.off,
    fill = "black"
  ), linetype = "dashed") +
  theme_bw() +
  labs(x = "TP53", col = "annotation")
## Warning: Ignoring unknown aesthetics: fill
beeswarm.plot

# sample.names <- intersect(rownames(drug_meta), colnames(exprMat))
# feature.df <- cbind(drug_meta[sample.names,], t(exprMat[,sample.names]))
# feature.df <- muvis::data_preproc(feature.df)
# feature.df <- feature.df[,-5183]
# mf <- muvis::min_forest(feature.df)

7 Differential Analysis with Limma

drug_nut <- diff_count(cps.data, drug_meta[which(drug_meta$TP53 == 0),], "nut")
drug_TP53 <- diff_count(cps.data, patMeta, "TP53")
meth_nut <- diff_count(meth.df, drug_meta[which(drug_meta$TP53 == 0),], "nut")
meth_TP53 <- diff_count(meth.df, patMeta, "TP53")
expr_nut <- diff_count(t(exprMat), drug_meta[which(drug_meta$TP53 == 0),], "nut")
expr_TP53 <- diff_count(t(exprMat), patMeta, "TP53")
top_sel(drug_nut)
## [1] "Nutlin.3a"   "RO5963"      "Cytarabine"  "Fludarabine" "Onalespib"  
## [6] "Palbociclib"
top_sel(drug_TP53)
## [1] "Nutlin.3a"   "RO5963"      "Fludarabine" "Cytarabine"  "Selinexor"  
## [6] "PF.3758309"
top_meth_nut <- top_sel(meth_nut, 0.05)
top_meth_TP53 <- top_sel(meth_TP53, 0.05)
top_expr_nut <- top_sel(expr_nut, 0.05)
top_expr_TP53 <- top_sel(expr_TP53, 0.05)
top_drug_nut <- top_sel(drug_nut, 0.2)
top_drug_TP53 <- top_sel(drug_TP53, 0.2)



vennList <- Venn(list(methylation_nut = top_meth_nut , methylation_TP53 = top_meth_TP53, expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))

upset(fromList(list(methylation_nut = top_meth_nut , methylation_TP53 = top_meth_TP53, expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53)), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")

8 Drugs

drug.annotation <- read.csv("../data/TP53/targetAnnotation_all.csv", sep = ";")
drug.annotation$nameEMBL2016 <- gsub(x = drug.annotation$nameEMBL2016, pattern = " +", replacement = ".")


drug.annotation_nut <- drug.annotation[drug.annotation$nameCPS1000 %in% top_drug_nut | drug.annotation$nameEMBL2016 %in% top_drug_nut,]
targets_nut <- strsplit(as.character(drug.annotation_nut$target), ", ")
targets_nut <- unlist(targets_nut)

drug.annotation_TP53 <- drug.annotation[drug.annotation$nameCPS1000 %in% top_drug_TP53 | drug.annotation$nameEMBL2016 %in% top_drug_TP53,]
targets_TP53 <- strsplit(as.character(drug.annotation_TP53$target), ", ")
targets_TP53 <- unlist(targets_TP53)


vennList <- Venn(list(drug_nut = targets_nut , drug_TP53 = targets_TP53))
Vennerable::plot(vennList, doWeights = F, show = list(Faces = FALSE))

upset(fromList(list(expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53, drug_nut = targets_nut, drug_TP53 = targets_TP53)), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")

genes <- Reduce(union, list(expression_nut = top_expr_nut, expression_TP53 = top_expr_TP53, top_meth_nut, top_meth_TP53))

# sample.names <- intersect(rownames(drug_meta), colnames(exprMat))
# feature.df <- cbind(drug_meta[sample.names,], t(exprMat[genes,sample.names]))
# feature.df <- muvis::data_preproc(feature.df)
# muvis::min_forest(feature.df[,-c(823:dim(feature.df)[2])])

9 Gene Set Enrichment Analysis

gmts <- list(HALLMARK=system.file("extdata","h.all.v5.1.symbols.gmt",package="BloodCancerMultiOmics2017"),
             ONCOGENIC=system.file("extdata","c6.all.v5.1.symbols.gmt", package="BloodCancerMultiOmics2017"),
             KEGG=system.file("extdata","c2.cp.kegg.v5.1.symbols.gmt", package = "BloodCancerMultiOmics2017"))

DEs <- list()
DEs$expression$nut <- expr_nut
DEs$expression$TP53 <- expr_TP53
DEs$methylation$nut <- meth_nut
DEs$methylation$TP53 <- meth_TP53
pathways <- DEs %>% map(function(x) x %>% map(function(y) names(gmts) %>% map(function(z) enricher(y, z, z))))

pathways <- names(DEs) %>% map(function(x) names(DEs[[x]]) %>% map(function(y) names(gmts) %>% map(function(z) enricher(DEs[[x]][[y]], paste(x, y), z)))) 
pathways[[1]][[1]] %>% map(function(x) plot(x[["enBar"]]))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
pathways[[1]][[2]] %>% map(function(x) plot(x[["enBar"]]))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
pathways[[2]][[1]] %>% map(function(x) plot(x[["enBar"]]))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
pathways[[2]][[2]] %>% map(function(x) plot(x[["enBar"]]))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
resss <- 1:2 %>% map(function(x) 1:2 %>% map(function(y) 1:3 %>% map(function(z) path.sel(pathways[[x]][[y]][[z]]$rg[[1]]))))
ress <- resss %>% map(function(x) x %>% map(function(y) unlist(y)))


vennList <- Venn(list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]]))
Vennerable::plot(vennList, doWeights = F, type = "ellipses", show = list(Faces = FALSE))

upset(fromList(list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]])), order.by = "freq", number.angles = 30, point.size = 3.5, line.size = 2, sets.bar.color = "blue")

Reduce(intersect, list(methylation_nut = ress[[2]][[1]] , methylation_TP53 = ress[[2]][[2]], expression_nut = ress[[1]][[1]], expression_TP53 = ress[[1]][[2]]))
## [1] "HALLMARK_APOPTOSIS"              "HALLMARK_UV_RESPONSE_DN"        
## [3] "HALLMARK_ALLOGRAFT_REJECTION"    "MTOR_UP.N4.V1_DN"               
## [5] "SIRNA_EIF4GI_UP"                 "KEGG_HEDGEHOG_SIGNALING_PATHWAY"